%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----  Simulate the market allocation under a centralizing platform  -----%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Setup workspace
close all
clear 

machine = '/Users/giulia/';

set(0,'defaultTextInterpreter','latex')

M = 15;

g_euler = 0.57721;

scale = 10^5;

%% Load matching function estimates

load(path_mf)
options = optimset('Algorithm','interior-point','display','on','MaxIter',90000,'MaxFunEvals',90000,'TolFun',10^-9,'TolX',10^-9,'TolCon',10^-9);

cobb_par = zeros(15,2);

for i =1:15
    col = freights(:,i)~=0;
    s = ships(col,i);
    f = freights(col,i);
    m = matches(col,i);
    lb = [0,0];
    ub = [Inf,0.95];
    theta0 = [rand(1),rand(1)];
    [cobb_par(i,:),~,~] = fmincon(@(theta)cobb(theta,f,s,m),theta0,[],[],[],[],lb,ub,[],options);
end

s_values = [1,linspace(0.95,1.95,10)];
v_values = linspace(0.75,5,11);

clearvars -except s_values v_values scale machine cobb_par

%% Load cost estimates

load(path_ships);

load(path_exports);

%% Platform Allocation
% We next examine whether a centralizing platform can bring us closer to the first-best. 
% By centralizing transactions, these platforms facilitate meetings and can improve market efficiency. 
% However, a centralizing platform does not act like a benevolent planner. Instead, it chooses prices to maximize its profit.

%
% Initialize the starting point
%

s_i = s_eff;
q_ij = m_eff.*Gdest_eff;
q_i = nansum(q_ij,2);
radius = 0.0;
delta_0 = (1 + unifrnd(-radius,radius,15,1)).*(q_i./s_i).*nansum(Gdest_eff.*meanprices_eff,2);
Gdest_0 = q_ij./nansum(q_ij,2);
lambda_0 = nansum(q_ij,2)./s_i;

%
% Setup convergence parameters
%

alpha = 10^-5;
coeff_q = 10^-6;
coeff_s = 10^-6;
iter = 0;
W_save = 0;
Niter = 10000000;

%
% We iterate over the matrix q_ij of carries traveling loaded on each route
% (ij) and over the number s_i of carriers waiting empty
%
while (iter<Niter)&&(diff>exp(-3))
    
    %
    % given a guess for q_ij and s_i, we compute the destination matrix G, 
    % matching rate lambda, and matches q_i
    %

    Gdest_0 = q_ij./nansum(q_ij,2);
    lambda_0 = nansum(q_ij,2)./s_i;
    q_i = nansum(q_ij,2);
    
    %
    % compute the average payoff delta_1 and value functions that implement
    % q_ij and s_i
    %

    [U_k,W_k,J_k,delta_1,inner_iter] = inner_loop(s_i,q_i,Gdest_0,lambda_0,delta_0,alpha,ksi,cs,cu,sigma,g_euler,S_obs,M,beta);
    
    % if the algorithm to compute delta_1 doesn't converge, choose a
    % smaller smoothing parameter

    if inner_iter == 10000
        alpha = alpha./10;
        [U_k,W_k,J_k,delta_1,inner_iter] = inner_loop(s_i,q_i,Gdest_0,lambda_0,delta_0,alpha,ksi,cs,cu,sigma,g_euler,S_obs,M,beta);
        [alpha, inner_iter] %#ok<NOPTS>
    end
    
    %
    % compute the gradient of the monopolist's profits with resepct to q_ij 
    % and s_i
    %
    
    [dpi_dq,~,~] = Dprofit_Dq(q_ij,beta,delta,s_i,cobb_par,E,M,r,K_obs,c_inv,W_k,U_k,J_k,sigma,sf);
    [dpi_ds,~,~] = Dprofit_Ds(q_ij,beta,delta,s_i,cobb_par,E,M,r,K_obs,c_inv,cu,sigma,W_k,U_k,J_k,sf);

    %
    % updte the guesses
    %
    coeff_mat = coeff_q.*ones(size(q_ij));
    coeff_mat(q_ij<0.01) = 10^-8;
    q0_ij = q_ij +coeff_mat.*dpi_dq;
    s0_i = s_i +coeff_s.*dpi_ds;
    
    diff = max(abs([dpi_dq(:)./(1 +q_ij(:));dpi_ds(:)]));
    
    iter = iter+1;
    q_ij = q0_ij;
    s_i = s0_i;
    delta_0 = delta_1;

end

%
% Save the allocation
%

q_m = nansum(q_ij,2);
Gdest_m = q_ij./nansum(q_ij,2);
ccp_m = ccp_equil(W_k,M);
s_m = s_i;
lambda_m = nansum(q_ij,2)./s_m;
f_m = inverse_m(q_m,s_m,cobb_par,M,sf);
lambda_f_m = nansum(q_ij,2)./f_m;
etilde_ij = entrants(q_ij,s_i,cobb_par,M,delta,sf);
Gunc_m = zeros(size(Gunc_ineff));
Gunc_m(:,2:16) = etilde_ij./E;
Gunc_m(:,1) = 1-nansum(Gunc_m,2);
J_m = J_k;
U_m = U_k;
W_m = W_k;
delta_m = delta_1;
meanprices_m = prices(q_ij,beta,delta,s_m,cobb_par,E,M,r,K_obs,c_inv,sf);
tot_ships_m = S_obs;
Jf_m = zeros(M,M);
profits = W1;
for i=1:M
    
    alpha = (1 - delta.*beta.*(1-lambda_f_m(i)));
    
    Jf_m(i,:) = lambda_f_m(i).*(r(i,:)-meanprices_m(i,:))-c_inv(i,:);
    
    Jf_m(i,:) = Jf_m(i,:)./alpha;
    
end
w_m = Welfare_new(Gdest_m,Gunc_m,ccp_m,s_m,f_m,q_m,cu./sigma,cs./sigma,c_inv,K_obs,ksi,r,E,sigma,Jf_m,W_m./sigma,U_m./sigma,sigma_exp,S_obs);
